Study of Vortex Avalanches using the Bassler-Paczuski Model 



G. Mohler and D. Stroud 
Department of Physics, The Ohio State University, Columbus, Ohio 43210 

(February 1, 2008) 



Abstract 

We have carried out a model calculation of the flux noise produced by vor- 
tex avalanches in a Type-II superconductor, using a simple kinetic model 
proposed by Bassler and Paczuski. Over a broad range of frequencies, we 
find that the flux noise S^{uj) has a power-law dependence on frequency, 
S^{uj) ~ with s ~ 1.4 in reasonable agreement with experiment. In 

addition, for small lattices, the calculated has a high-frequency knee, 

which is seen in some experiments, and which is due to the finite lattice size. 
Deviations between calculation and experiment are attributed mostly to un- 
certainties in the measured critical current densities and pinning strengths of 
the experimental samples. 

I. INTRODUCTION 

Bassler and Paczuski Q have recently proposed a model which simulates 
the development of the Bean critical state in a type-II superconductor, as 
well as of other types of so-called self-organized critical behavior In such 
a system, vortices are injected into one edge of an initially empty sample, 
and are free to fall off the opposite edge only. As vortices are driven into the 
system from the "loading edge," they push already-present vortices further 
into the sample. These existing vortices then pile up in such a way that the 
vortex density gradient approaches a critical value which is determined by 
system parameters such as pinning strength and density. Once the gradient 
reaches this critical value, the system is said to have achieved a "self-organized 
critical state." 

Of particular interest in this context, is the development of "avalanches" 
from the Bean critical state. When the flux density gradient exceeds its critical 
value, the injection of even a single excess vortex into the system can start 
a chain reaction of vortex motion, known as an avalanche, which may have 
a very large scale in both space and time and which causes the gradient to 
relax back to its critical value. Vortex avalanches are typically characterized 
by their duration, their linear extent in the direction of average vortex motion, 
and the number of vortices forced from their original positions ( "topplings" ) 
by the avalanche. These and other avalanche characteristics are expected to 
obey various scaling laws Q. 
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In earlier numerical models, the generation of vortex avalanches was typi- 
cally studied numerically at a "microscopic" scale, i. e., at a length scale where 
individual vortex displacements were calculated from certain postulated force 
laws In practice, this type of model imposes severe constraints on the size 
of system which can be studied numerically; such constraints in turn make 
it difficult to study the critical behavior of vortex avalanches. But since an 
avalanching system is expected eventually to achieve a self-organized critical 
state, the large-scale behavior of the system should be describable without 
tracking the motion of individual vortices. The BP model takes advantage of 
this expectation by focusing attention only on the large-length scale behavior 
of this system, which can be modeled using only a coarse-grained lattice. This 
choice makes feasible the simulation of large systems. 

II. THE BASSLER-PACZUSKI MODEL: DEFINITION AND METHOD OF 

CALCULATION 

In the Bassler-Paczuski (BP) model, one considers a distribution of vor- 
tices on a two-dimensional simple hexagonal lattice (see Fig. The scale of 
the lattice is assumed to be such that the vortices can be treated as point-like 
objects which exist entirely on the lattice sites. Each lattice site is capable 
of containing multiple vortices, and the number of vortices per lattice site is 
a kind of coarse-grained vortex density |^ . Each lattice site is also described 
by a pinning potential, which is chosen at random from a suitable distribu- 
tion. Periodic boundary conditions are imposed at the top and bottom of 
the lattice. Vortices are injected into the left-hand edge, and removed from 
the right-hand edge, as described in more detail below. As more and more 
vortices are injected into the lattice, the repulsion between vortices eventually 
overcomes the attraction exerted on the vortices by the pinning sites. As a 
result, the vortices slowly migrate from left to right across the lattice. When 
they reach the right-hand edge, in this model, they are assumed to "fall off" 
the lattice and are removed from the system. Vortices are forbidden to exit 
from the left side of the system, or to re-enter the lattice from the right-hand 
side once they fall off. 

An important parameter in the BP model is the rate of vortex injection. 
In the "slow-driving limit" IQ, a vortex is injected into the lattice only after 
an avalanche has ceased. This slow-driving limit is similar to that studied in 
models which treat the microscopic (i. e., the individual vortex) degrees of 
freedom together with explicit force equations for each vortex. For example, 
Nori et al who have studied such a model, introduce a new vortex into 
the lattice only when the lattice-averaged vortex velocity has fallen below a 
certain threshold Q. In the present paper, we inject vortices either slowly or 
quickly, so that we may make predictions for either a hypothetical slow-driving 
regime or the high magnetic field ramping rates studied in most experiments 
i- 



2 



In the BP model, the force acting on a vortex at site x in the direction of 
site y is taken to be 

Fx^y = —V{x) + V{y) + [m{x) — m{y) — 1] + r[m{xl) + m{x2) — m{yl) — m{y2)]. (1) 

Here xl and x2 are the two sites (other than site y) which neighbor x on the 
hexagonal lattice; yl and y2 are the two sites other than x which are nearest 
neighbors to y; V{z) is the strength of the pinning potential on site z; and 
m{z) is the vortex population at site z (taken to be a positive integer), r is 
the strength of the repulsion between vortices on neighboring sites, in units 
such that the on-site repulsion strength is normalized to one. 

The BP model is defined by three parameters, denoted r, p, and q. r is the 
nearest-neighbor repulsion parameter mentioned above, which is assumed to 
be smaller than unity; p is the pinning strength; and q is the probability that 
a particular site contains a nonzero pinning potential. Thus, any particular 
site has a pinning strength of p or with probability q and 1 — q respectively. 
The distribution of pinning sites is determined at random once for the entire 
lattice, at the start of the simulation, and is held constant thereafter. 

A vortex move is carried out as follows. First, the force on a vortex 
is calculated for each of the three possible directions of vortex motion. If 
exactly one of the three forces is greater than zero (i. e., acts towards a 
nearest neighbor site), then one vortex is moved by one lattice spacing in 
the direction of that force. If more than one force is greater than zero, the 
direction of motion is chosen randomly from the set of positive forces. This 
calculation is carried out for each site containing at least one vortex, and all 
the vortex populations are updated in parallel - that is, each site is examined 
once in a given time step, and at most one vortex is moved from a particular 
site during a time step. 

To carry out the calculations, we constructed a simple hexagonal lattice 
using the unit cell shown in Fig. |l], with cell length and cell height Uy = 
\^ax- The total lattice size was chosen so as to have dimensions LxUx x Lytty, 
as shown in the Figure, with Lx = ^Ly, and Ly integer. For this choice of 
dimensions, and for Lx > 32, the rectangular sample was found to be wide 
enough to prevent individual avalanche events from overlapping one another 
in space. For a given lattice size, we carried out simulations of avalanche 
formation and evolution for a range of pinning strengths and vortex injection 
rates. 

III. RESULTS 

A. Numerical Results Obtained from BP Model 

All of the runs were carried out using 2^^ ~ 10^ time steps. For most 
choices of parameters and vortex injection rates, it was found that the vortex 
population quickly reached a steady state of the form predicted by the Bean 
model That is, the magnetic induction gradient dB/dx from left to right 
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(i. e., in the direction of vortex injection) rapidly approached a constant value 
(see Table). According to the Bean model, this constant value is related 
to the critical current density Jc of the superconductor (see below). 

Our goal in this work is to calculate the flux noise produced by vortices 
falling off the right edge of the sample. This goal is similar to that of Jensen, 
who has also used a coarse-grained model to study the power spectrum of a 
self-organized critical state [^]. Jensen's model, however, differs from the BP 
model in several respects. For example, the Jensen model forbids occupancy 
of a site by more than one vortex, and it also allows vortices to fall off both 
the right and left edges of the system. Thus, we expect the two models to 
give different predictions for noise spectra. 

To calculate the flux noise in the present model, we used the following 
approach. In experiment, the flux noise is measured by placing a detector coil 
in the center of a hollow sample ring as shown in Fig. |3|. According to 
Faraday's Law, the voltage generated in the coil is proportional to the rate of 
change of vortex population within the coil area. To display this, we plot in 
the Figures the Fourier transform 

5$(/) = Limr^ool C ^'{t) exp(-i27r/i)dtp, (2) 
Jo 

where N'{t) represents the rate of change of vortex population in the region 
to the right of the sample at time t. (Our transverse periodic boundary 
conditions correspond to a ring geometry.) In our calculations, N'[t) is taken 
as the number of vortices which fall off the right-hand edge of the sample 
in a single time step at time t. We evaluate the Fourier transform using 
the integral given above, then smooth the resulting power spectrum using a 
Savitzky-Golay filter as described in Ref. |12|. We used first order smoothing 
over twenty points on either side of the data point. 

Our results are shown in Figs. ^ and ^. In general, for all the spectra 
shown, the power spectrum exhibits a power-law dependence on frequency, 
over a portion of the frequency range, i. e., S^{f) ~ resembling that seen 
in the experimental results of Field aZ In Fig. s ~ 1.0 for spectra a, 
b, and c. s ~ 1.4 for d. In general, s is found to increase from 1.0 to 1.4 when 
the injection rate is reduced below one vortex per time step across the entire 
left-hand edge. However s never achieves the experimentally observed value 
of ~ 1.5, as seen, for example, in the results of Q. 

For the relatively small lattices, L < 32, we often see a characteristic fea- 
ture in the power spectrum which is missing at larger sizes. Namely, for an 
injection rate of one vortex per time step, where the vortex can enter any- 
where along the left-hand edge, the power spectrum exhibits a plateau; see 
Fig. ^(a). This plateau disappears as the injection rate is decreased. For 
example, at an injection rate of 1 vortex per 100 time steps [cf. Fig. |5|(6)], no 
plateau is apparent. The plateau is due to the emergence of a new length scale 
in the lattice for small lattices and large injection rates, namely, the small- 
est linear dimension of the lattice itself. Under such conditions, the lattice 
size limits the avalanche dimensions. This new length scale in turn produces 
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the high-frequency "knee" in the noise spectrum. Specificahy, at these fre- 
quencies, lattice-wide avalanches dominate, drowning out the high-frequency 
noise produced by sporadic independent avalanches. As in the other power 
spectra generated by the model, the calculated slope b = —dln{Sy{(jj)) / dln{uj) 
for frequencies to above the knee of Fig. ^ is smaller than the experimental 
value hexp', our calculated value in this regime is 5 ~ — 1 compared to the 
experimental value of b 

exp ^ — ^ ' 

A vortex injection rate greater than 1/1 (i. e., greater than one vortex per 
lattice site per time step) was also attempted, so that multiple vortices were 
randomly injected into the lattice simultaneously. The goal of this test was to 
check if a "knee" could be created at larger length scales at sufficiently high 
injection rates. But this test resulted in a power spectrum with nearly white 
noise, i. e., 6 ~ 0. It is possible that, with the rather limited lattice size of 
L = 32, the size of avalanches is inevitably restricted purely by size limita- 
tions, whereas in the larger lattices, the vortices exhibit a different avalanche 
pattern, unconstrained by lattice size. 



B. Connection to Experiment 

In order to compare our model results with the experiments of Field et al, 
we need to estimate the lattice constant and time step of the simulation. To 
obtain the lattice constant, we look at a cross-section of the lattice, as seen 
in Fig. 1^. Then integrating Ampere's Law, 

^ ^ 47r , 
V X B = —J 

c 

along the path in Fig. |2[ and assuming a constant current density equal to 
the critical current density, J = Jc (as expected in the Bean critical state), 
we find 

If we multiply this expression by the area of a narrow vertical strip of the 
lattice, jUyUx, we obtain the flux through the vertical strip centered at x and 
of width 



$(x) = m - (^a,a,) 



To proceed further, we write where Ux is the distance between two 

opposite sides of the hexagonal cell (see Fig. Q), and we denote by m{x) the 
total flux through a column parallel to y, of width a^, and centered at Uxa, 
Then it follows from eq. |5| that 



X- 



minx) = m(0) - J 



Finally, if the cells are equilateral hexagons of side a, then Ux = a\/3 and 



Qy = 3a, and hence 

f ^ fn\ ( ^^JcLa^ \ 
m{nx) = m(0) - I — — I 

We can use this equation to interpret our numerical results for the vortex 
densities and flux noise spectrum. We also need an estimate of Jc, the critical 
current density for the sample studied experimentally. In this case, the sample 
is the composite material NboArTio^^ (typically in the form of a solid solution 
with precipitates which constitute the pinning centers). Experiments suggest 
a value in the range of 3 x 10^^ statamperes/cm'^ ]ll|, from which we can 
estimate a from 

m(0)$oc y/^ 
OttLVc ) 

Some representative values of a, as calculated from our simulations, are shown 
in Table I. 

Once a is determined, the simulation time constant, tq, can be estimated 
by equating the Lorentz force per unit length on a vortex to the pinning force 
per unit length, i. e. 

fpin — fb 

We assume that fpin can be interpreted as a vortex drag force proportional to 
the average vortex velocity i. e., fpin = rjv, where r] is an effective friction 
coefficient. Then writing out the Lorentz force on a unit length of a single 
vortex explicitly gives: 

rjv = 



Now a vortex can move only one lattice constant per time step, so that v ~ 
a/rp. Furthermore, r/ = where B is the local magnetic induction, p is 

pc 

the flux flow resistivity, and c is the speed of light. Combining these relations 
gives 



Ba 

To 



pcJc 

If we assume a magnetic induction of 5 kG, a lattice constant of 3.22 x 10~^ cm 
as suggested by the estimates in the Table, a critical current of 3 x 10^^ 
statamps/cm^ as is typical for the experiments of Field et al, and a flux flow 
resistivity of 1.11 x 10^^^ sec, as suggested by experiments on NbTi alloys 



|11|, then one obtains a time constant tq ~ 2 x 10 ^ s. We now discuss the 



implications of these estimates. 
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IV. DISCUSSION 



If the above value of tq is substituted into our calculated frequency spec- 
trum, we find that the spectrum has a power-law frequency dependence over 
a slightly different frequency range than that in which power-law behavior 
is seen by Field et al in Nbo ^jTio^^^ (cf. Figs. ^ and ||). Furthermore, the 
strength of the model flux noise power is smaller than the experimental val- 
ues by a factor ranging from two (at high frequencies, for high injection rate) 
to five (at low frequencies, for low injection rate) orders of magnitude. 

These differences between our calculated power spectrum and the mea- 
sured spectrum are not surprising. The model is a only a rough approximation 
to a real material, since it uses an artificial kinetics rather than a more realis- 
tic (and more computer-intensive) dynamics for the calculation. Furthermore, 
because of the finite size of the simulation samples, they produce significantly 
less flux noise than would be generated by a real material. Nevertheless, our 
model does give the qualitatively correct behavior: it leads to a power-law 
exponent which seems to approach the observed value for a sufficiently low 
vortex injection rate and a sufficiently large lattice. 

The noise spectra deviate from power-law behavior at low frequencies. 
Specifically, they all show a weak peak at low frequencies, followed by a fur- 
ther decrease with diminishing frequency and finally an increase at still lower 
frequencies. We believe that this peak occurs near the frequency of a char- 
acteristic "lattice resonance." Qualitatively, this frequency is the ratio of a 
characteristic length, i. e., the linear dimension of the lattice, and a char- 
acteristic time, which is the time required for the lattice to "reload" to its 
Bean critical state between avalanches. For a lattice of size Lx = 64, having 
pinning parameters (0.1, 5.0, 0.1) and an injection rate of one vortex per time 
step, the resonant frequency occurs at about 3Hz (using our estimates for 
lattice constant and time step). For the same parameters but slower injection 
rate (one vortex per IOOtq), the noise spectrum still shows a peak but now 
at a lower resonance frequency of 0.3 Hz. For other values of the pinning 
parameters and injection rates, the resonance becomes less conspicuous. The 
resonant peak also becomes much less conspicuous in a larger lattice. In this 
case, the lattice simply contains more vortices. Since some of these may not 
join the primary avalanche, they tend to move separately, reducing the ap- 
parent reload time and hence producing a stronger low-frequency noise signal. 
The same argument holds when the pinning strength is increased at constant 
lattice size. In this case, there are once again more vortices in the system 
than at weaker pinning, and hence, a stronger low-frequency flux noise signal. 
Finally, at high injection rates (e. g., in the L^. = 32 lattice at an injection 
rate of one vortex per time step), the vortex system is in a state of continual 
avalanche motion, since this rate is very large for such a small lattice. As a 
result, the power spectrum, which reflects that of the avalanches, is relatively 
flat as a function of frequency. 

Finally, we comment briefly on the relation of our calculation to a similar 
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study done by Nori and collaborators ||T^. In contrast to our work, their 
calculations were carried out using dynamical equations for individual vor- 
tices (assuming overdamped motion and a particular force law to describe 
vortex- vortex interactions; the equations were solved using molecular dynam- 
ics methods. Their calculations were carried out solely in the slow-driving 
limit; the resulting noise spectra, like those presented here in the slow-driving 
limit, exhibited a power law in agreement with experiments carried out in that 
slow-driving regime. But our calculations do have the advantage of speed: a 
realistic MD model using overdamped dynamics of individual vortices requires 
approximately 10^ hours on a machine with parallel processing. By contrast, 
our power spectra were calculated using a single Digital Alphastation 255. 
for only 10^ hours. This speed allows us to vary the parameters extensively, 
in particular examining the effects of different driving rates, which in turn 
permits investigation of the shoulders in the power spectra mentioned above. 
Of course, the kinetic approach does have the compensating disadvantage of 
invoking a drastic simplification to the true equations of motion; but it still 
appears to preserve much of the relevant physics. 

To summarize, we have shown that a relatively simple kinetic model of 
vortex avalanches near the Bean critical state [Q] gives rise to flux noise which 
qualitatively resembles experiment, without using computer- intensive dynam- 
ical equations for individual vortices. Specifically, the model predicts a power 
law flux noise spectrum with approximately the correct power-law exponents, 
as well as a length-scale-induced knee in the power spectrum for a sufficiently 
small lattice at high frequencies. Although the model calculations differ in 
detail from experimental results, they do show many qualitative similarities 
and trends in the context of rather simple kinetic equations. 
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Table I: Calculated lattice constant a and magnetization density m(0) for several 
lattice sizes and two choices of model parameters 



L 


m(0)($o) 


a(10 ^cm) 


slope^ 


128" 


7500 


3.22 


1.83 


64" 


1900 


3.24 


1.86 


64^ 


2300 


(3.45) 


2.25 


32" 


510 


3.37 


1.99 



"parameters (r, p, q) = (0.1, 5.0, 0.1) as defined in text; 
''parameters (0.1, 12.0, 0.1). The lattice constant for b was calculated using 
the same critical current as for a. In actuality, the larger pinning strength 
would call for a larger critical current. 

^the slope is calculated from the columnar vortex density versus distance into 
the sample, in terms of unit cells; i.e. {m{Lx)/Ly — m{0)/Ly)/{Lx — 0). 
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FIGURES 



FIG. 1. Hexagonal lattice, with a rectangular overlay displaying the unit cells used in the 
calculations. Our "unit cell" is a rectangle enclosing four lattice points as shown. 

FIG. 2. Cross-section of the vortex population in the sample (schematically indicated by the 
length of the vertical arrows) after the vortex density has relaxed to the steady-state vortex profile, 
showing the linear decrease in magnetic induction with distance into the sample, as predicted by 
the Bean model. The path of the line integral used to integrate Ampere's Law is shown in boldface. 

FIG. 3. Schematic of experimental arrangement for noise measurements, as viewed from above. 
The arrows denote the average direction of vortex motion, under the influence of an increasing 
external magnetic field applied to the outside of the torus. The coil is used to measure the flux 
noise. 

FIG. 4. Calculated voltage noise spectra due to flux motion into the interior of the torus, for 
several choices of parameters and primarily for larger lattices, r and q are 0.1 for all spectra. The 
legend indicates the pinning strength, p, and the lattice width (in units of Ox-.) The injection rate is 
one vortex injected per time step (fast injection) for spectra a, b, and c, and one vortex injected per 
100 time steps (slow injection)for spectrum d. The numerical data were collected for a run of 2^" 
time steps, using a Savitzky-Golay smoothing filter of order 1, using 40 points. The translation into 
real frequency is made assuming the estimates of time constant and lattice constant as described 
in the text. 

FIG. 5. Additional calculated vortex noise spectra, primarily for smaller lattices. These results 
show a possible "knee" in curve a. This curve represents an injection rate of one vortex per time 
step, which would correspond to experimental data taken at a high rate of vortex injection. This 
curve should be contrasted to curves b and c, which are taken at lower rates of magnetic field 
ramping (1 vortex per 100 times steps) and show no such knee. 
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Power Spectra, Larger Lattices 
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Power Spectra, Smaller Lattices 

r=0.1 ,p=5.0,q=0.1 for all spectra 
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